ON THE REGULARIZATION OF FREDHOLM INTEGRAL 
EQUATIONS OF THE FIRST KIND* 
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Abstract. In this paper the problem of recovering a regularized solution of the Fredholm integral 
equations of the first kind with Hermitian and square-integrable kernels, and with data corrupted by 
additive noise, is considered. Instead of using a variational regularization of Tikhonov type, based 
on a priori global bounds, we propose a method of truncation of eigcnfunction expansions that can 
be proved to converge asymptotically, in the sense of the L 2 -norm, in the limit of noise vanishing. 
Here we extend the probabilistic counterpart of this procedure by constructing a probabilistically 
regularized solution without assuming any structure of order on the sequence of the Fourier co- 
efficients of the data. This probabilistic approach allows us to use the statistical tools proper of 
time-series analysis, and in this way we attain a new regularizing algorithm, which is illustrated 
by some numerical examples. Finally, a comparison with solutions obtained by the means of the 
variational regularization exhibits how some intrinsic limits of the variational-based techniques can 
be overcome. 
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1. Introduction. We consider the Fredholm integral equations of the first kind 

(1.1) (Af)(x)= I" K(x,y)f(y)dy = g(x) (a<x<b) 

J a 

whose kernel K(x J y) is supposed to be Hermitian and square integrable; i.e., 

(1.2) K(x,y) = K(^x) 
and 




(1.3) / < / \K(x,y)\'dx }dy< oo 



Then A : L 2 (a, b) — > L 2 (a, b) is a self-adjoint compact operator. 

For simplicity we shall suppose hereafter that the kernel K, the function g, and the 
unknown function / are real- valued functions; in addition, we assume that the interval 
[a, b] is a bounded and closed subset of the real line. 

The Hilbert-Schmidt Theorem guarantees that the integral operator A admits a 
set of eigenfunctions {tpk}T an d, accordingly, a countably infinite set of eigenvalues 
{^fe}i°- The eigenfunctions form an orthonormal basis of the orthogonal complement 
of the null space of the operator A and therefore an orthonormal basis of L 2 (a,b) 
when A is injective. For the sake of simplicity only this case will be considered, 
although this assumption can be easily relaxed with slight technical modifications. 
The Hilbert-Schmidt theorem also guarantees that lim; £ _, 00 Xk — 0. Furthermore, we 
shall suppose hereafter that the eigenvalues are ordered as follows: \\ > A2 > A3 > .... 
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In view of the Hilbert-Schmidt theorem we associate with the integral equation (|l.lfl 
the following eigenfunction expansion: 

where gt = (g,ipk), ((v) denoting the scalar product in L 2 (a,b)). The series (11.411 
converges in the sense of L 2 . 

Remark. If the support of the data does not coincide with that of the solutions, 
i.e., A : L 2 (a, b) — ► L 2 (c, d) with [a, b] different from [c, d), the problem can be worked 
out in terms of singular values and singular functions of the operator A , and all of 
the following results can be easily reformulated. 

In view of the fact that there always exists some inherent noise in the data, instead 
of l|l.lf> we have to deal with the following equation: 

(1.5) Af + n = g (g = g + n), 

where n represents the noise. Therefore, instead of expansion (|1.4(l we have to consider 
the following expansion: 

(1.6) 

where g k = (g, ij)k). Expansion 1)1.6(1 is generally diverging because ~g does not belong, 
in general, to the range of the operator A. This is precisely a manifestation of the 
ill-posed character of the Fredholm integral equation of the first kind. 

Several methods of regularization have been proposed (see |10l 1141 HB] and refer- 
ences therein); all of them modify one of the elements of the triplet {A, X, Y}, where 
A is the integral operator defined by 1(1. 1(1 . whereas X and Y are, respectively, the 
solution and the data space (in our case X = Y = L 2 (a.b)). Among these methods 
the procedure, which is probably the most popular, consists in admitting only those 
solutions that belong to a compact subset of the solution space X. In particular the 
famous method of Tikhonov leads to the construction of "regularizing operators" by 
the minimization of "smoothing functionals" . In this latter functional the smoothing 
term is obtained precisely by restricting the admitted solutions to a compact subset 
of the space X: then the continuity of A -1 follows from compactness. This restriction 
is realized by the use of a priori bounds which can be written assuming some prior 
knowledge of the solution. Therefore, in addition to the inequality 

(1-7) \\Af-g\\<e 

which corresponds to a bound on the noise (|j • || denoting the norm in L 2 (a,b)), one 
also considers an a priori bound on the solution of the following form: 

(1-8) \\Cf\\ z <E, 

where Z denotes the "constraint space" and, accordingly, C is the "constraint opera- 
tor" . From the bounds ((1.7(1 and 1(1.8(1 we are led to define the regularized solution as 
the minimum of the following functional: 

(1.9) $(/) = \\Af- gf + « 2 \\C.f\\l , (a = (-|)) . 
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In spite of several significant merits, this procedure is not free from defects. Con- 
cerning the possibility of writing suitable a priori bounds on the solution, we want to 
remark strongly that two different types of problems must be distinguished: 

a) synthesis problems; 

b) inverse problems, 

and to note that both are frequently solved by the use of Fredholm integral equations 
of the first kind. In the first class of problems, that basically consists in finding the 
source that produces a prescribed effect (e.g. prescribed boundary values), the a priori 
bounds are intrinsic of the problem itself, whereas this is not always the case for the 
second class. As typical examples we can consider: 

a') the antenna synthesis; 

b') the signal recovery. 
The problem of the antenna synthesis consists in determining, within a certain degree 
of approximation, the current intensity that generates a desired radiation pattern. It 
can be formulated in terms of Fredholm equation of the first kind |181 01) and, conse- 
quently, it presents the typical pathology of the ill-posed problems. In this problem 
the a priori bound on the ohmic losses associated with the current intensity is neces- 
sary and can be regarded as a natural constraint intrinsic of the problem. Conversely, 
in the case of the signal recovery problem, the a priori bounds can be written only 
if prior knowledge on the signal is assumed. Generally, it is possible to have some 
a priori information regarding, for instance, the support of the signal or requiring 
the function representing the signal to be nonnegative. But even in these cases the 
prior knowledge could be insufficiently specific to be peculiar of the function to be 
reconstructed, and arbitrary, though reasonable, constraints must be added to solve 
the problem. Strictly connected with this question there is the crux of the matter: 
the practical choice of the regularization parameter a (see formula for a fixed 

g, when the a priori bound (|1.8(l is unknown or it is not sufficiently precise. 
Moreover, let us note that the functional 1)1.9(1 works as a filter whose action is smooth- 
ing the Fourier components g k for high values of k. But it is easy to exhibit examples 
of signals whose Fourier components are small, or even zero, for low values of k, while 
the significant contributions of the signal are brought by those components at inter- 
mediate values of k, which are smoothed out by the action of the filter. In these 
situations the standard regularization method fails, showing that the only existence 
of the minumum of functional (|1.9fl does not guarantee the bulk of the signal had 
been really recovered. This delicate point will be illustrated with numerical examples 
in section 

We suggest a different approach which is based on the following observation: 
for the moment, suppose that the moduli of the noiseless Fourier coefficients \gk\ 
are monotonically decreasing as k increases; then, although the formal series l|1.6(l 
diverges, nevertheless the effect of the error remains limited in the beginning of the 
expansion, and there exists a point (a certain value of k) where divergence sets in. 
Thus, the idea is to stop the expansion at the point where it turns to diverge. This 
rough and qualitative description can be put in rigorous form by proving that even 
if the series (|1.6|l diverges, nevertheless it converges (in the sense of i 2 -norm) as 
e (i.e., the bound on the noise) tends to zero. This result, which has been proved 
by two of us (see JZj), does not give (except in very particular cases) a practical 
numerical method for finding out the truncation point (i.e., the value of k) where to 
stop expansion However, here we prove a probabilistic generalization of the 

results presented in |17| by removing the quite restrictive assumption that the Fourier 



4 



E. DE MICHELI, N. MAGNOLI, AND G. A. VIANO 



coefficients |<7fc| of the signal to be recovered are monotonically decreasing. Compared 
to J7| the significance of the new results is relevant. First, the hypothesis made 
in ^7] on the order of the coefficients |<7fe| leads to a regularization procedure that 
essentially works as an ideal low-pass filter, and, as previously discussed, this does 
not guarantee to recover correctly the signals whose bulk is localized at intermediate 
frequencies. Conversely, in this paper it will be shown how to construct a regularized 
solution without assuming any kind of order on the coefficients \gt\ by exploiting the 
tools supplied by the information theory. This result will lead to a more effective 
regularizing algorithm which is based on a suitable statistical analysis of the data and 
whose main feature is indeed the frequency selectivity. Second, from the application 
point of view, the hypothesis on the order of the coefficients \g k \ is too restrictive; 
thus, by removing it, a much larger class of real signals can be practically analyzed. 
These questions are precisely the contents of sections [21 and 0] We will prove, indeed, 
in Section[3]that it is possible to split the noisy Fourier coefficients g k into two classes: 

i) the Fourier coefficients g k from which a significant amount of information on 
fk = (f,ipk) can be extracted; 

ii) the Fourier coefficients g k that can be regarded as random numbers because 
the noise prevails on the coefficients gk- 

In section 0] it will be shown how it is possible to separate practically the co- 
efficients g~ k into these two classes by the use of statistical tools supplied by the so 
called "time-series" analysis. Therefore, we can practically construct an approxima- 
tion which converges to the real solution, and furthermore we can have some confidence 
that the bulk of the function / has been effectively recovered. 

The paper is organized as follows. In the first part of section [3 a short sketch of 
the variational method based on the minimization of functional i|1.9|) is given. This 
will be done in order to have explicitly the formulae which will be used in section 0] 
where our procedure and the variational one will be compared. The second part of 
section is devoted to the probabilistic formulation of the regularization problem in 
a quite general setting. In section [21 we start illustrating the asymptotic convergence 
of the eigenfunction expansion (in the sense of L 2 -norm) as e tends to zero; then this 
result is reconsidered from the viewpoint of probability and information theory. Here 
a key role will be played by the Bayes formula: it will provide the various terms of 
our approximation, which will be proved to be a probabilistically regularized solution 
of ijl.lfl . The first part of section 01 is devoted to the discussion of the statistical 
tools that are necessary for practically recovering the regularized solution from finite 
samples of noisy data. Finally, some numerical examples are given in the second part 
of section 0] 

2. Variational and probabilistic regularization. 

2.1. Variational regularization. After the classical book of Tikhonov and Ars- 
enine [2'6\ the literature on the theory and applications of the variational regularization 
has been rapidly growing (see, for instance, |14j1. In order to compare our algorithm 
with this classical one, some formulae and results of the variational regularization will 
be here recalled here (see El GDI G2J for proofs and details). 

Let us characterize, first of all, the constraint operator C and, accordingly, the 
constraint space Z. Let us take a constraint operator C such that C*C and A* A com- 
mute (this assumption does not restrict the theory and the applications significantly 
[HIE]). Then, the space Z is composed by those functions / G L 2 (a,b) such that 
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||C/||z is finite; i.e., 

OO 

(2.1) \\Cf\\% = (C*Cf, f) = 5>£|A.| 2 < ex.. 

k=l 

Now we consider the ball Uz = {/ € Z | J2kLi c k\f k \ 2 — an d * ne restriction Ao 
of the operator A (see Ql-lfl ) to the ball Uz- Then, the following propositions can be 
proved. 

Proposition 2.1. //linifc^ooC 2 = +oo the operator Aq 1 is continuous. 
Proposition 2.2. The functional $(/), with a — (e/E), has a unique minimum 
which is given by: 

(2.2) /, = [A*A+(-|) 2 C*C]- 1 A*g. 



By expanding g in terms of ifck (eigenfunctions of the operator A), we have: 

oo . 

k=i A k + c fc ("eJ 

Next, we have the following proposition. 

Proposition 2.3. The following limit holds true for any function f satisfying 
the bounds and 

(2.4) lim 11/ - A|| =0 (E fixed). 



In numerical computations it is often convenient to use truncated approximations. 
For instance, one can derive from the smoothed solution (|2.3|) the following truncated 
approximation: 



(2.5) 



fc=i 



where k a is the largest integer such that 
(2.6) X k > (|) 



pa- 



Proposition 2.4. The following limit holds true for any function f satisfying 
bounds and 



(2.7) 



lim 

e— 



/ - ft 



(1) 



= (E fixed). 



In several problems a weaker a priori bound should be used by setting C = I (the 
identity operator). Therefore, instead of bound l|1.8|) . we have 

/ oo \ V2 

(2.8) ll/H =(ElM 2 J ^ K 
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In this case the unique minimum of functional (|1.9|) is given by 



(2.9) 



/i 2) - E 



t4\ 



2 f/tl 



k=i K + (§) 

and, accordingly, the following truncated approximation can be introduced: 



(2.10) 



fe=i 



A3) 9k i 



where A;^ is the largest integer such that 



(2.11) 



A fc > — . 
_ E 



Both fi 2 ^ and /| 3 ' ) converge to / as e — > in a weak sense. In fact, as shown in 
[23 IS] ! the following proposition can be proved. 

Proposition 2.5. For any function f which satisfies the bounds \1.7j and ^2. St . 
the following limits hold true: 



(2.12) 



lim 



([/-/i 2) ] ,v)\ =0 (||«|| < 1, E fixed), 



(2.13) 



lim 



([/-/i 3) ] ,v)\ = (||«|| < 1, .E /ized). 



2.2. Probabilistic regularization. Here we want to reconsider ()1.5[) from a 
probabilistic point of view. With this in mind we rewrite (|1.5fl in the following form: 

(2.14) A£ + C = tj, 

where £, C and f?, which correspond to /, n and g respectively, are Gaussian weak 
random variables (w.r.v.) in the Hilbert space L 2 (a,b) [2]. A Gaussian w.r.v. is 
uniquely defined by its mean element and its covariance operator; in the present case 
we denote by R^, Rcc and R^ the covariance operators of £, £ and r\ respectively. 
Next, we make the following assumptions: 

I) £ and C have zero mean; i.e. = = 0; 
II) £ and ( are uncorrelated, i.e. R^^ = 0; 
III) R7^ exists. 

The third assumption is the mathematical formulation of the fact that all the com- 
ponents of the data function are affected by noise. As it is shown by Franklin (see 
formula (3.11) of JJ|), if the signal and the noise satisfy assumptions I) and II), then 

(2.15) R m = AR^A* + R cc 
and the cross-covariance operator is given by 

(2.16) = RtfA*. 

We also assume that R^^ will depend on a parameter e that tends to zero when the 
noise vanishes; i.e., 



(2.17) 



R C C = <?N, 
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where N is a given operator (e.g., N = / for the white noise). 
Now we are faced with the following problem. 

Problem. Given a value g of the w.r.v. r\ hnd an estimate of the w.r.v. £. 

A linear estimate of £ will be any w.r.v. £^ = Lrj, where L : Y — » X, is an 
arbitrary linear continuous operator. Then, from a value g oi rj one obtains the linear 
estimate of the w.r.v. £. Now a measure of the reliability of the estimator L is 
given by 

(2.18) 5 2 (e 1 v;L) = E{\(t-L V ,v)\ 2 }, (v e X = L 2 (a, &)), 

where E{-} denotes the expectation value. Then, we have the following proposition. 

Proposition 2.6. If the covariance operator has a bounded inverse, then 
there exists a unique operator Lq that minimizes S 2 (e,v;L) for any v € X, and it is 
given by 

(2.19) L = R iv R-j = %A* [AR^A* + R^]' 1 . 



Proof See HE]. □ 

The w.r.v. LqT] is called the best linear estimate of £, and, given a value g of 77, 
the best linear estimate fi 4 ^ for the value of £ is 

(2.20) /i 4) = 9, (A* = A). 

v ' J AR^A* + i?cc 

If ^ and Lr/ have finite variance, then the global mean-square error may be defined as 
follows: 

(2.21) S 2 (e,L)^E{U-Lr 1 \\ 2 }. 

When the operator L which minimizes (|2.18|l does exist, it also minimizes the global 
error (|2.21l) if Lqt\ has finite variance; i.e., if Tr (LqR^Lq) < 00, then the following 
proposition can be proved. 

Proposition 2.7. If the following assumptions 

i) is an operator of trace class; 

ii) Rqq = e 2 N has bounded inverse; 

iii) the equation Af — 0, where f € Range ^i?^ 2 ^ , has only the trivial solution 
/ = 

are satisfied, then the following limit holds true: 

(2.22) lim(5 2 (e) = 0, 

where S 2 (e) = inf L S 2 (e; L) . 
Proof. See HEj. □ 

Let us note that S 2 (e) = S 2 (e; Lq) when Lq does exist and is unique. 

If we want to compare the probabilistic results obtained above with the variational 
ones, which have been obtained by the use of eigenfunction expansions, we must 
expand £ and ( in terms of the eigenfunctions of the operator A (i.e. {^fe}? )- Their 
Fourier components are the random variables = (£,ipk) an d Cfc = (CiV'fc), whose 
variances are given respectively by pf. and e 2 v\. Next, in addition to the assumptions 
I)TII) made before, we make the following hypothesis in spite of the fact that it turns 
out to be completely unrealistic (see section 0}: 
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IV) the Fourier components of £ are mutually uncorrelated as well as the Fourier 
components of (. 

Therefore, if is bounded (i.e. sup^l/e 2 ;/ 2 ,) < oo), then the operator Lq exists 
and the best linear estimate (|2.20(l can be written as 

(2-23) /i 4) = £ AV % ai 

k=l k ^ k k 
Finally, the quantities S 2 (e,v;Lo) and 5 2 (e) become 
6 2 (e,v;L ) = E{\(£- L r),v)\ 2 } = 



(2.24) 



00 22 

Pk v t ,„, ,2 



= (1% - LoRr,r,L^]v,v) =e 2 J2 \2 n i\\l„i \ Vk 

fc=1 A kPk + 6 v k 

and 

°° n 2 u 2 

(2.25) <5 2 (e) = £ 2 (e; L ) = Tr [i? 55 - L i^i$] = e 2 £ ' * 



\2 „2 , 2 ..2 ' 
fc=1 A fePfc + 6 ^fe 



and we have the following proposition. 

Proposition 2.8. The following statements hold true: 

i) for any v G X (X = L 2 (a, b) ) 

(2.26) limS 2 (e,v;L ) = 0, 

£— >0 

ii) ?/ TrR^ < 00, then 

(2.27) lim<5 2 (e) = 0. 

£^0 



3. Information theory and regularization. 

3.1. Asymptotic convergence, in the L 2 -norm, of the eigenfunction ex- 
pansion. In the variational regularization, use is made of global a priori bounds (e.g., 
formulae (|1.8|) or (|2.8[l l. which are the natural constraints in the case of synthesis prob- 
lems where the variational approach is certainly appropriate. But these global bounds 
are not necessarily given in the case of inverse problems where the prior knowledge 
on the solution can be, in several cases, rather poor. Moreover, in the truncated solu- 
tions derived by the methods of variational regularization, the point at which to stop 
the expansion is obtained by comparing the eigenvalues with the ratio (e/E) (i.e., 
formula ^.Hfl ). or with (e/E)\ck\ (see formula (|2.6|0 . In both cases this approach 
appears quite unnatural from the viewpoint of the experimental or physical sciences, 
whose methodology would rather suggest to stop the expansions at the value ko of k 
such that for k > ko the Fourier coefficients g^ of the noiseless data are smaller or 
at most of the same order of magnitude of e, and, consequently, it is impossible to 
extract information from the corresponding coefficients ~g k . With this in mind, and 
assuming that the noise is represented by a bounded and integrable function n(x) 
which satisfies the following condition: 



(3.1) 



sup \n(x)\ < e, x £ [a,b], 
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the following results have been proved by two of us: 
Lemma 3.1. The following statements hold true: 



(3.2) =ll/ll 2 = C i (Ci = constant), 



k=l 



oo /_ s 2 

I 9k 



(3.3) 2^ lit' = ify^RangeiA), 



k=l 



(3.4) \\mg k =g k , Vfc. 

£—►0 



// ko (e) is defined by 



(3.5) feo(e) = max <j m e IN : ^ ( |M < d j> , 



m /_ \ 2 



fc=l 



f/je 



(3.6) lim fco(e) = +oo 



Proof. See [Hj. □ 

Now we can introduce the following approximation 

fco(e) _ 

(3.7) / (£) = E g H k 

k=l k 

and prove the following theorem. 

Theorem 3.2. The following equality holds true: 



(3.8) Urn f-fi 



f (e) 



= 0. 



Proof See [I7j. □ 

If we consider a sequence of noisy data g which tends to g for e — > in the sense 
of the L 2 -norm (i.e., lim e ^o \\g — g\\ = 0), then f^ will tend to / as e — > in the sense 
of the L 2 -norm (i.e., lim e ^ ||^ e) -/|| = 0). In fact, since ||5-g|| 2 = J2T=i l<?fc-5fc| 2 , 
the lim^o ||<? — g\\ — implies that for any k, lim e ^o5fc = 9k, and in view of Lemma 
13. II and Theorem 13. 21 it can be concluded that lim e ^o ll/d^ ~ /II = 0. Therefore, from 
approximation i|3.7[) we can derive an operator B defined by: 

fcoO) _ 

(3-9) Bg=J2^ k , 

fc=i k 

which continuously maps (i.e., preserving the convergence) the data ~g into the solution 
space X . Thus, continuity has been restored without requiring compactness. 
Two types of difficulties still remain: 
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a) how to determine numerically the truncation point fc (e), if the norm of the 
function / (i.e., the constant C\ — ||/|| 2 ) is unknown; 

b) in any case the convergence of approximation l|3.7|) is not sufficient to guar- 
antee that the bulk of the unknown function / has been really recovered. 

We can give a satisfactory answer to these questions only in very specific and peculiar 
cases, as we will explain below. Suppose that the moduli of the Fourier coefficients |<jffe| 
are monotonically decreasing for increasing values of k. Since g k = gk + n k , it turns 
out that at a certain value ko of k we have \gk\ — \n k \ < e. The Fourier coefficients of 
the noiseless data are of the same order of magnitude as the Fourier components of 
the noise, and at this point we cannot extract any information from the noisy Fourier 
coefficients ~g k . Let us now introduce the function M(m) = JZfcLi (Vk/^k) 2 , whose 
relevant properties are: 

1) It is an increasing function of m; 

2) If e is sufficiently small and the values of \g k \ are monotonically decreasing 
for increasing fc, M(m) presents a "plateau" when it reaches the value C\. 
Indeed, from formula .(j|> in Lemma l3.1l it follows that M(m) remains nearly 
constant when it attains the value C\. An explicit numerical example of this 
"plateau" is given in Figure ETTt ) in section 

This "plateau" corresponds to the order-disorder transition in the coefficients g k : for 
k < /co(e) the data g k prevail on n k whereas for k > fco(e) the noise components n k 
are larger or, at least, of the same order of magnitude of the noiseless data. However, 
it must be remarked that in practical cases to single out the plateau which does really 
correspond to the order-disorder transition in the coefficients g k can be made difficult 
by the presence of other spurious plateaux due to the erratic behavior of the noise. 
Furthermore, if the coefficients g k are negligible for low values of fc, and the actual 
bulk of information is located only at intermediate values of fc, there could be no 
numerical evidence of such a plateau in spite of the fact that the convergence guaran- 
teed by Theorem 13. 21 remains true. Then we are forced to look for other methods that 
overcome these difficulties. This issue will be investigated by means of probabilistic 
methods, as will be illustrated in the next subsection. 

3.2. Bayes formula, information theory, and regularization. Here our 
goal is to find a probabilistic extension of the result of Theorem 13.21 in which the 
assumption requiring the Fourier coefficients \g k \ to be monotonically decreasing will 
be removed. In fact, we will show how to construct a regularizing solution from the 
noisy data, disregarding the order of the coefficients \gk\- For this purpose, we turn 
(|2.14|) into an infinite sequence of one-dimensional equations by means of orthogonal 
projections: 

(3.10) Aft£ fc + 0fc=f?fc, (k = 1,2,...), 

where = (£,ipk)> Ck — (C^fc), i] k = (r],ip k ) are Gaussian random variables. Here 
we retain assumptions I)-III) made in section |2~2*1 but we remove assumption IV). 
In fact, there is no reason to assume that the basis {^k]T which diagonalizes the 
operator A also diagonalizes the covariance operators R^, R^(, R m US]. Therefore, 
we can introduce the variances p\ = (R^ipk,ipk), £ 2 v k — {Rcc^kt^k), tfcPk + = 
(R^rfipk, "0/c), without assuming that the Fourier components of £ (and analogously 
also for ^ and rj k ) are mutually uncorrelated. In view of the assumptions I) and III) 
the following probability densities for £ k and Cfc can be assumed: 

(3.11) Ml) = -^« p {-(£!)}, (,= 1,2,...) 
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and 

By the use of the (|3.1U[) we can also introduce the conditional probability density 
p Vk (y\x) of the random variable r\ k for fixed !; k = x, which reads 

Pm{y\ x ) = -75= — ex P 



r^tvk { 2e 2 



'k 



(3.13) = _^ eX p{-^^ 



Now let us apply the Bayes formula that provides the conditional probability density 
of £fc given rj k through the following expression: 

(3.14) ViM v) = Va f^ X) - 

Thus, if a realization of the random variable r\ k is given by g k (see the formulation of 
the problem in section IT2"|) . formula l|3.14|l becomes 

(3.15) p ik (x\g k ) = A fc exp|-^||exp|-^-| - |^ \ (A k = const.). 

Now the amount of information on the variable which is contained in the variable 
r\k can be evaluated. We have ^Sj 

(3.16) = -ilog(l-r2), 
where 

(317) r 2_ |E{&^.}| 2 (AfePfe) 2 



Thus, 



(3-18) m,Vk) = l log fl + 44 

From equality Ij3.18|l it follows that J(£k,i]k) < ^log 2 , if AfcPfc < ev k . Thus, we are 
naturally led to introduce the following sets: 

(3.19) 3 k = {k : X k pk > e^fe} , 



(3.20) N k = {k : X kPk <eu k }. 

Reverting to the conditional probability density (|3.15|) . it can be regarded as the 
product of two Gaussian probability densities: p\ (x) — A k exp {-x 2 /2p 2 k } and 
p 2 (x) = A[ 2) exp {-(A 2 /2e 2 ^) ( x _ (g fe /A fc )) 2 }, (A k = A k 1] ■ A{ 2) ), whose variances 
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are respectively given by pk and (eVk/^k)- Let us note that if k € 3^, the variance 
associated with the density p 2 (x) is smaller than the corresponding variance of pi (x) , 
and vice versa if k E Therefore, it is reasonable to consider as an acceptable ap- 
proximation of the mean value given by the density P2{x) if k € Jfe, or the mean 
value given by the density p\ (x) if k G We can write the following approximation: 

(3.21) (&)=J ft {ke3kh 

{ o (fceNfc). 

Consequently, given the value ~g of the w.r.v. t], we are led to consider the following 
estimate of £: 



(3.22) 




However, these are only heuristic considerations based on plausible arguments. They 
will become rigorous statements only if it will proved that they lead to a solution B~g 
which is probabilistically regularized. For this purpose, the global mean-square error 
associated with the operator B, i.e., E — £??7|| 2 j, must be evaluated, and we have 
the following proposition. 
Proposition 3.3. 

i) 7/limfc_ >00 (AfcPfc/z/ / t) = 0, then the set 3k is finite for any fixed positive value 
oft; 

ii) assuming that the limit stated in i) holds true, and, in addition, that is 
an operator of trace class, then the following relationship holds: 

(3.23) E{U-B V r}= £ ^+£^#<°°- 

fceN* fc£J fc k 



Proof. The proof of statement i) is obvious if we recall the definition of the set 3k 
(formula Ij3.19|l h Statement ii) follows easily from the equality 

(3.24) E — Br]f\ = Tr (R a - R^A*B* - BAR a + BR m B*) 

and by the use of formulae f2~l"5|) . (|2~T7jl . and l|3.22[l . □ 

In order to prove that approximation (|3.22(l is regularized, we need the following 
auxiliary lemma. 

Lemma 3.4. Let fcy(e) be defined as follows: 



(3.25) 




where T = TrR^ is finite. Then the following statements hold true: 
(3.26) i) limfc 7 (e) = +00, 

! k -< 2„ 2 °° 1 

E pi =°- 
fe=i * fc=fe~+i 



REGULARIZATION OF FREDHOLM INTEGRAL EQUATIONS 



13 



Proof, i) Let fc 71 denote the sum {k 1 + 1). Then suppose that the limit 
does not hold. This latter assumption would imply that there exists a finite number 
M, which does not depend on e, such that fc 7l < M, Furthermore, this bound 
should remain true for any sequence e, tending to zero. Then, we have the following 
inequalities: 



Now for any sequence Ci tending to zero, we have 

M oo 

(3.29) r<]T p 2 <Y j pI = t, 

k=l k=l 

which is contradictory. Then, limit (|3.26() holds. 

ii) From Ylk=i Pk ~ Tri?^ = T < oo, and in view of statement i), it follows that 
lim e ^ J2kLk., Pk — 0- Regarding the sum X)fc=i( e " ''^k I '^D » we can P rocee d as follows. 
From formula Q3.25|l we have 

k k 

(3.30) E^# + E^<r. 

fc=l fe k=l 

Then 

(3.3D E^ r -E^= E ^ 



a/ 



but in view of the fact that lim e ^o J2'kLk^ 1 Pk = 0, we have lim c _>o Sfcli( e2jy fcMl) = 
0, and statement ii) is proved. □ 

We can now prove the following theorem. 

Theorem 3.5. If the covariance operator is of trace class, and if the set 3 k 
is finite (see Provosition \3.3\l . then the following limit holds true: 

(3.32) limd 2 (e,B) = lim - Br)\\ 2 } = 0; 



i.e., approximation \3. 2ty is probabilistically regularized. 

Proof. In view of formula H3.23fl in Proposition 13.31 the proof of equality (|3.32l) 
reduces to the proof of the following limit: 

(3.33) ^{E^+E^W 

UeJ fc k k€N k ) 

Regarding the first sum of (|3.33J) . we divide the set 3k into two subsets defined by 

(3.34) 5^ = {k G 3 k : k < fc 7 }, 

(3.35) 3 k 2) = {k G 3 k : k > k 7 }; (3 k = 3 ( k 1] U 3 ( k 2) ). 
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Accordingly, we can write 



(3.36) E^=E^+E 



\A \2 1 \2 

^ Afc ke^ Afe ke^ h 



Then J^, -a) ( e2j/ fcM 2 .) — zCfe=i( e2 ^fc/^fc)' anc ^ m y i ew of Lemma IXH (where we 

zee j fc 

proved that lim e ^o IZfeLi ( e2 ^l/^I) = 0) it follows that 

2 2 

(3.37) lim V 1-^=0. 

ke3 



(i) - k 



Regarding the term J2k^i i2> ( e2 ^fcMfc)> smce k &3k then p 2 > (e^/A 2 ) and therefore 



2,, 2 



3C 



2 



(3-38) E ^ E 



\2 

,(2) A fe 



But, as we have seen in Lemma f3. 41 lim £ ^o /C^Lfc Pfc = 0> an d consequently 

2 2 

(3.39) lim V 1-^=0. 

c^o ^— ' A 2 



<2) - fc 



We can conclude that lim e _ Sfeea fc ( e2 ^fcMfc) = 0- Regarding the sum X^fceN p\, we 
proceed in an analogous way by splitting the set Jik m to two subsets defined by 

(3.40) = {k G K fe : fc < fcy}, 

(3.41) 5 = {^^:fc> fc 7 }; (W fc = U W< 2) ). 
Accordingly, we write 

(3.42) E ^ = E ^ + E p2- 

If fc e N^ 3 and by the use of inequality p\ < (e 2 ^ 2 /A 2 ) (because k € N&) we can write 



(3-43) J2 P t<j2 



A 2 



But in Lemma we proved that lim e -_> YlkLii^^k/^k) = 0> an< ^ therefore we have 
lim e ^ St, C T\r (1) Pk = 0- Regarding the second term on the right-hand side of formula 
we have 



(3-44) E ^ E & 



hex* 

„2 



But, again, lim e ^ J2T=k 11 Pk = °> an d tnen nm e^o L^in' 2 ' Pfe = °- E 
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Remarks, i) It is worth it to notice that the proof of Theorem 13.51 does not re- 
quire any type of order in the sum i|3.22[l . In fact, the only assumption that is a 
strictly decreasing sequence does not evidently imply that the terms (X k p k /ev k ) have 
any type of monotonicity in k, and, consequently, the sum 13.22fl cannot, in general, 
be regarded as an ordered sum of terms up to a certain maximum value of k. Thus, 
unlike the regularized solutions l|2.3jl . (|2.5(l . 12.9(1 . Ij2.10jl and also (|3.9|) . Bg features 
frequency selectivity, which is obtained by evaluating the information content of the 
noisy Fourier coefficients. 

ii) Notice that the estimate (|3.22|) associated with the operator B represents a proba- 
bilistically regularized solution, in the sense of the formula l|3.32[l . even if, in general, 
it does not minimize the global mean-square error (|2.21|) . 

At this point in order to apply the results of this section, statistical methods that 
allow for splitting the coefficients ~g k into the sets Ok and N/c must be investigated. 
These methods will be illustrated in the next section. 

4. Numerical analysis: the regularizing algorithm. 

4.1. The correlation function of the noisy data. The application of the 
results of the previous section to a Fredholm equation of the first kind would involve 
using statistical tools for the determination of the two sets 3 k and Nfc. In this section 
this issue is discussed and the basic steps of a numerical algorithm for constructing 
the regularized solution Bg from the noisy data g are outlined. For simplicity we 
shall work throughout only with data corrupted by white noise. However, provided 
the independence assumption between £ and £, more general cases involving "col- 
ored" noise could be treated by using suitable methods, for instance, "prewhitening" 
transformations 0, whose discussion is beyond the scope of this section. Here our 
goal is to show that statistical estimates of the amount of information carried by the 
Fourier coefficients ~g k can be sufficient to construct a satisfactory regularized solution. 
Furthermore, the direct comparison of the numerical results clearly evidentiates how 
some inherent limitations of the variational regularization scheme are overcome. 

Following the analysis of the previous section, we are now faced with the problem 
of separating the Fourier coefficients ~g k into two classes; one containing all the Fourier 
coefficients of the noisy data which are correlated, the other containing the ~g k that 
can be regarded as random numbers. This task can be achieved by computing the 
correlation function of the random variables rj k : i.e., the probabilistic counterpart of 
the coefficients g k : 

an A (u E{[m, -E{y kl }}[Vk 2 -E{ Vk2 }}} 

[ ' ^U*2) E{[vki _ E{mA] 2 V /2 E{[vk2 _ E{vk2WV/ 2- 

In practice, just a finite realization {g k }^ of the random variables r] k is available, 
from which estimates &g of the autocorrelations can be obtained by regarding the data 
{Hk\i as a finite length record of a stationary random normal series. In principle the 
assumption of stationarity of the series {rj k } is not correct because in general the 
moments of the random variables r] k will depend on fc, but from the practical point 
of view this is usually the only possible chance. In fact, in many areas of application, 
it is difficult or even impossible to have multiple independent realizations {g^}^ of 
the process {r/fc}, so estimates of ensemble averages cannot be computed. Thus, we 
are forced to introduce the working hypothesis that the process {r) k } is stationary 
in wide sense [5], that is A r) (fci,/c2) = A r/ (ki — £2), and compute the estimates of 



16 



E. DE MICHELI, N. MAGNOLI, AND G. A. VIANO 



the autocorrelation coefficients by means of the ergodic relation between ensemble 
and time (i.e., the index k in our case) averages. Of course, such a restriction can 
be removed whenever many independent sets of data {g^.} 1 ^ would be available for 
evaluating ensemble averages. Anyway, we will see later in the discussion of the 
algorithm how an ambiguity in the reconstruction of the regularized solution Bg due 
to the assumed invariance for fc-translation of {rjk} will be removed. 
A number of estimators of the autocorrelation function have been suggested by statis- 
ticians and their properties are discussed in detail in |15| . An estimate which is widely 
used by statisticians, and in the following examples as well, is given by 

N-n 



E (5k ~ (9k))(9k+n - (9k+n)) 



(4.2) %(n) = ±± — , n = 0,...,N-l, 



N—n N—r 



1/2 : 



k=l k=l 



where 



N—n N— r 



( 4 - 3 ) (9k) = Jf—^ E 9k, (9k+n) = Jf—^ E 9k+n- 

k=l k=l 

Equation 14. 2[) . which is based on the scatter diagram of ~g~k+n against g k for k = 
1,..,N — n, represents the maximum likelihood estimate of the autocorrelation co- 
efficients of two random variables r/k and r\k+n whose joint probability distribution 
function is bivariate normal. 

In order to identify the structure of the series {gk}\ so that we can separate 
correlated components from the random ones, it is necessary to have a crude test on 
whether <%(n) is effectively zero. It has been shown by Anderson that the distri- 
bution of an estimated autocorrelation coefficient, whose theoretical value is zero, is 
approximately normal. Thus, on the hypothesis that the theoretical autocorrelation 
A,j(ri) = 0, the estimate Sj(n) divided by its standard error asiji) will be approx- 
imately distributed as a unit normal deviate. This fact may be used to provide a 
rough guide as to whether theoretical autocorrelations are essentially zero. To this 
purpose it is usually sufficient to remember that, for normal distribution, deviations 
exceeding two standard errors in either direction have a probability of about 0.05, so 
that the 95% confidence interval of the estimate is approximately Sj(n) ± 1.9617,5(71). 

Estimated autocorrelations can have rather large variances and can be highly 
correlated with each other J3J I12j , so that care is required in the interpretation of in- 
dividual autocorrelations. In particular, moderately large estimated autocorrelations 
can occur after the theoretical autocorrelation function has damped out and, in any 
case, it must be considered that an estimated autocorrelation function always exhibits 
less damping than the theoretical one, as the estimated autocorrelations are inflated 
by sampling fluctuations (see also the following Example 1). Thus, in order to avoid 
a purely empirical analysis of the autocorrelations, it is necessary to assume a rough 
model of the series that allows to evaluate the order of magnitude of the sampling 
errors as(n) associated to the autocorrelation estimator. 

According to the discussion of Section l3~2l since we are expected to find the set 3k 
to be finite, we are also expected that the autocorrelation function A j; (n) will vanish 
beyond a certain lag tiq. Thus, in what follows, it will be assumed that there exists 
an index no such that A j; (n) = for n > Uq. In this case, if the record length N is 
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large enough (i.e., such that 0(1/N 2 ) terms can be neglected), use can be made of the 
Bartlett's approximate expression for the variance of the estimated autocorrelations 
of a stationary normal process |3] : 

1 f n ° 1 

(4.4) var[%(n)] ~ |l + 2^ A^(«) j , for n > n . 

To use l|4.4fl in practice, the estimated autocorrelations &g are substituted for the 
theoretical ones A, ( , and in this case we shall refer to the square root of (|4.4I) as the 
large-lag standard error cr${n; no) \7\. 

The index no is actually recovered in a recursive way through an hypothesis 
generation-verification procedure. Starting from the assumption that the series is 
completely random, i.e., no — 0, the standard error ag(n; 0) is computed and the first 
index n > such that |<%(n)| > 1.96 175(71; 0) is searched for. If there exists such 
an index n, it becomes the new candidate to be n , i.e., we set no = n, as(n;no) is 
computed, and again it is tested whether the series is compatible with the hypothesis 
that A n (n) — for n > uq. The whole procedure is repeated until no new index n is 
found. Formally, no is then defined as 

(4.5) n = max{n > : Vn E {n 7 N- 1], | %(n) |< 1.96 a s (n,n)}. 

The set Q of the lags corresponding to autocorrelation values that are effectively 
different from zero and, consequently, indicating lack of randomness of the coefficients 
g k , is defined as: 

(4.6) Q = {0 < n < n :\ <%(n) |> 1.96 a s (n, 0)}. 

Let N c be the number of elements of Q. 

As previously discussed, as a consequence of the inevitable assumption of sta- 
tionarity of the process the Fourier coefficients ~g k that are correlated cannot be 

determined in a unique way from the set Q. In fact, an integer 6 Q just indicates 
a strong correlation between at least two Fourier coefficients ni apart. This means 
that, in principle, any couple (fffc^S^+nJ for any integer 1 < ki < (N — n^) could 
have generated such a strong correlation at the lag n^. Thus, from the set Q we can 
construct N c families Fi defined as 

(4-7) ^ = {(5 fel ,5 fcl+ „,)}if = T ) ' i=^-,N c 

from which the couples of coefficients ~g k that are likely to be correlated can be selected. 
In theory, that is for N — * 00, the N c indices ki and the N c elements n, e Q are 
mutually dependent. In fact, any two coefficients ~g ka i~g kfs which are selected from the 
families Fi must satisfy the pairwise compatibility conditions requiring \k a — kp\ S Q. 
Or, in other words, it can be seen that, given the set Q, the number Nj of admissible 
Fourier coefficients ~g k is combinatorially constrained to be 

(4.8) i (1 + y/l + 8N C ) < Nj < N c + 1. 

The left inequality in (|4.8|l follows directly from the observation that the maximum 
number of correlations among Nj coefficients is ( ^ ) , then N c < ( ^ ) , whereas the 
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right inequality expresses that at least (Nj — 1) distinct correlations can be computed 
among Nj coefficients (i.e. N c > Nj — 1). For instance, if N c — 2, we have from 
inequalities (|4.8|l that there need to be Nj — 3 coefficients g k to construct the set Q, 
or, referring to l|4.7|l . that the two indices k\ and k2 must coincide, i.e., k% = ki > 1. 
In any case, the compatibility conditions are not sufficient to constraint in a unique 
way the selection of the coefficients g k and, consequently the construction of the reg- 
ularized solution. 

In practice, that is when the record length N is finite and particularly when the signal- 
to-noise ratio (SNR) of the data ~g is small, the compatibility constraints cannot be 
assumed to be satisfied. In fact, because of the sampling fluctuations in the estimates 
5g(n), some correlations which are actually different from zero could be incorrectly 
detected by the procedure discussed above. However, we shall see later in the discus- 
sion of the numerical examples how the compatibility constraints can provide us with 
a confidence check on the reliability of the regularized solution Bg. 

In order to recover in a unique way from the set Q the Fourier coefficients that 
are likely to be correlated, we adopt the following criterion suggested by the definition 
itself of the autocorrelation function: for any m £ Q, i = 1, N c , we select the pair 
(~9k*i9k*+n ) giving the maximum contribution to the autocorrelation estimate Sjirii); 
i.e., we define k* as 

(4.9) fc* = arg max {\g k g k+ni \}, i = l,...,N c , 

and, accordingly, we can define the set of frequencies 3 k exhibiting correlated Fourier 
coefficients as 

(4.10) 3 k = {k*}? c U{k* + n i }?% 

where each element of 3 k is counted only once. 

4.2. Numerical examples. Throughout this section we shall consider as a sam- 
ple problem the integral equation with kernel 



< 1, 

< 1 



tsI ^ f 0-~x)V if 0<l/<I< 

(4-11) K(x,y) = < 'K ., ~ ~ ~ 

\ x (1 - y) if < x < y < 

whose eigenfunctions and eigenvalues are, respectively, 

(4.12) ip k (x) = \/2 sin(A:7ra;), 

(4-13) A fc = ^. 

The data g(x) have been noised by adding white noise n{x), simulated by computer 
generated random numbers uniformly distributed in the interval [— e, e] (see also |22j 
for a very preliminary numerical analysis of this problem). The examples shown here- 
after differ for the choice of the input signal f(x) and for the values of the noise bound- 
ary e, whereas the performances of the algorithm are evaluated by direct comparison 
of the reconstructed signal with the true signal f(x). In every example reported 
here, the approximations obtained through the variational scheme (see section 12.111 
are computed by setting the constraint operator C such that c k = k, (k = 1, 2, ...), 
the parameter e corresponding to the boundary on the noise equal to the dispersion 
of the noise D e (see (^3), and the boundary E on the solution equal to the norm of 
the unknown function, i.e., E = ||/(cc)|| (see (| 1 . HI) ^ . 



REGULARIZATION OF FREDHOLM INTEGRAL EQUATIONS 



19 



Kir 



,B 



f(x) 



M(m) r 



D 



1 1 . 1 1 1 1 , 1 

6 8 10 12 

n 



JJ L 

14 16 



Fig. 4.1. Example 1: f^x) = (1 - x) sin(3 sin(3a;)), t = 10" 4 , SNR ~ 25.7dB, N = 512. 
(A) Noiseless Fourier coefficients g^. (B) Modulus of the autocorrelation function. The horizontal 
dotted straight line indicates the 95% confidence limit 1.96 crg(n;0) for a purely random sequence. 
The solid curved line indicates the confidence limit 1.96 erg (n; 3), n > 3. From the analysis of Sg(n) 
we have Q = {1,2,3} and 3^ = {1,2,3,4}. (C) Regularized solutions. The solid line represents 
the actual solution fi(x). The dots represent the reconstruction B~g. The crosses represent the 
variational solution fj^ obtained by using c^ = k; k a = 8 (see equations 12.5V and 12. 6Xl ). (D) Plot 
of the function M(m) = ^2 k _ 1 (<?fc/^fc) 2 - Notice that the value of M(m) corresponding to its first 
plateau, i.e., approximately for 4 < m < 10, is about the squared norm of the true solution. 



In Figure PT~Tl the analysis of the sample function fi(x) = (1 — x) sin(3 sin(3x)) 
with noise boundary e = 1CP 4 is summarized. The global SNR, defined as the ratio 
of the mean power of the noiseless data to the noise variance, was SNR ~ 25.7dB. 
The function f\{x) is characterized by having the bulk of information localized in the 
first few values of k (see the related noiseless coefficients in Figure l4~TK ) so that 
we expect that also a variational solution could provide a satisfactory reconstruction 
of the input signal. Figure |4"TT 3 shows the behavior of the autocorrelation function 
&a{n) along with the two lines indicating the statistical confidence limits we used to 
discriminate whether the autocorrelations are essentially null. The dashed horizon- 
tal straight line represents the threshold that we would have under the hypothesis 
of purely random sequence {g^}, whereas the solid line represents the threshold cor- 
responding to the model of autocorrelation function of ideal damped type. In this 
example we found uq = 3, Q = {1,2,3}, and the autocorrelation at n = 4 was re- 
jected in spite of its quite large value (see formula l|4.fi|l ). The direct inspection of 
the values of evk and gk in repeated realizations showed that for k = 5 the noise was 
usually larger than the Fourier coefficient, confirming hence the result that the auto- 
correlation <%(4) was abnormally inflated by the large autocorrelations at n — 1, 2, 3. 
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According to the criteria (|4.9|) and l|4.10|l , the set of frequencies whose corresponding 
Fourier coefficients exhibit strong correlations is 3f. = {1,2,3,4}. It is worth noticing 
that in this case the elements of 3k satisfy all the compatibility constraints, i.e. any 
difference between elements of 3k belongs to Q, and Nj satisfies constraints (|4.8|) . 
This complete cross-consistency between Q and 3k gives a high level of confidence in 
the result of the whole analysis. In Figure WHO , the true function /i (solid line), the 
regularized solution Bg (crosses) and the regularized function /i 1 ^ (dots) are com- 
pared. The truncation point of obtained through the criterion (|2.6p . was a = 8. 
Figure 14. 1C shows how in this case both regularization methods lead to comparable 
results, which are quite satisfactory approximations of the "unknown" function f\. 
The plot of the function M(m), displayed in Figure l4~Tt ) confirms the correctness of 
the two approximations. In fact, it clearly exhibits a "plateau", ranging from about 
to = 3 to rn = 10, that corresponds to the order-disorder transition of the coefficients 
g~k- Then it could be argued that for any truncation point belonging to this "plateau" 
the truncated approximation will hold coefficients g k whose information content is not 
completely obscured by the noise. In every example discussed here, the regularized 
solutions /i 2 ^ and /i 3 ^ (see (|2.9ll and H2.10|l 1 have also been considered, providing in 
all cases worse results (not plotted). 

The second and third examples, shown in Figure EOl are quite simple but a little 
tricky, and show the deep differences between our approach and the variational one. 
They consist of a finite linear combination of, respectively, 3 and 10 basis functions 
ipk (see the legend for numerical details), and, indeed, they have been chosen as 
typical signals in which the bulk of the information is not grouped in a single block of 
consecutive low frequencies. In these cases, setting global constraints on the solution, 
such as in the variational methods, leads inevitably to a failure, which is clearly evident 
from Figure I4.2C .F. since the lack of selectivity necessarily causes the regularized 
solution /i 1 ' to contain pure noisy components. On the contrary, the selectivity 
achieved through the analysis of the autocorrelation function overcomes this limit. 
In both examples the analysis of the autocorrelation function (see Figure I4.2I AD) 
led to the correct selection of the components that carry information in spite of the 
quite small SNR (in the Example 2, SNR ~ 0.55dB). Referring to the Example 2 
depicted in Figure l4~2*K .B.C. it can be observed that all the compatibility constraints 
are indeed satisfied; however, it is worth to remark that, because of the sampling 
fluctuation of the estimates %(n), the autocorrelation <%(6) was not always detected 
in different realizations of the noisy data {g~k}i ■ I n these cases the set 3k, computed 
from the set Q = {4,10} missing n — 6, is still correct, i.e., 3k — {3,7,13}, even 
though one compatibility constraint is not fulfilled. 

A more complex example is shown in Figure l4~3l Following the trace of the pre- 
vious example, here we have the input function which is characterized by having 
the significant Fourier components grouped in different ranges of the k axis. Conse- 
quently, the Fourier coefficients gk that clearly emerge from the noise (in this example 
e = 10 -4 ) are quite sparse in the range 1 < fc < 12 (see Figure FOR 1 ). The plot of the 
regularized solution Bg, obtained from the analysis of the autocorrelation function 
shown in Figure FOB . shows an acceptable agreement with the real solution f^, even 
though the procedure failed in detecting the coefficient at k = 5. On the contrary, 
the "nontruncated" (in the sense that the sum runs up to N) solution /+ (see H2.3I) ). 
which is displayed in Figure l4~3"£ >. yields a rather poor reconstruction either because 
the constraint operator C smooths out too many frequencies or because distortions are 
introduced by those coefficients which are essentially noise (e.g., k = 2,6,7,8,9,10). 
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Fig. 4.2. Example 2: f 2 (x) = 5sin(37nr) + 10sin(77nc) + 15 sin(137rz), e = 310" 3 , 
SNR ~ 0.54AB, N = 512. (A) Modulus of the autocorrelation function. Q = {4,6,10}, 
3fe = {3,7,13}. (B) Comparison between the actual solution f2(x) (solid line) and the regular- 
ized solution B~g(x) (dots). (C) Comparison between the actual solution f2(x) (solid line) and the 
approximated solution (x) with k a = 9 (see criterion 12. 61 ) ). (D) Example 3: Modulus of the 
autocorrelation function. fs{x) = ^^—l a i s ' m (kjirx) , with aj = {17, 23, 27, 33, 43, 55, 68, 70, 77, 81} 

and kj = {5,9,13,17,18,23,24,25,31,33}. e = 10~ 3 , SNR ~ 9.79dB, N = 1024; Q = 
{4, 5, 8, 9, 12, 13, 14, 15, 16, 18, 19, 20, 22, 24, 26, 28}, 3 k = {5, 9, 13, 17, 18, 23, 24, 25, 31, 33}. (E) Com- 
parison between the actual solution f3(x) (solid line) and the regularized solution Bg(x) (dots). (F) 

Comparison between the actual solution $2{x) (solid line) and the approximated solution f^(x) 
with k a = 27. 



Of course, the variational reconstruction could be considerably improved by choosing 
a more appropriate operator C and different values for the parameters e and E, but 
this would require more precise a priori knowledge on the actual solution. 

In conclusion, some final remarks. The method of regularization based on the 
analysis of the correlation function of the data allows to pick out the Fourier com- 
ponents of the noisy data which are likely to carry exploitable information on the 
unknown solution, and at the same time, for rejecting the ones dominated by the 
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Fig. 4.3. Example 4: f 4 (x) = (1 - x) sin(5sin(12x)), e = 10" 4 , SNR ~ 4.6dB, N = 
512. (A) Noiseless Fourier coefficients g^. (B) Modulus of the autocorrelation function. Q = 
{1,2,3,5,7,8,9}, 3fe = {1,3,4,9,11,12}. (C) Comparison between the actual solution fn(x) (solid 
line) and the regularized solution Bg(x) (dots). (D) Comparison between the actual solution fi(x) 
(solid line) and the variational solution f*(x) (see I2.3tl ). 



noise. Frequency selectivity is not featured by methods of regularization that basi- 
cally work as low-pass filters, and we have seen this inherent limit through examples 
in which frequency selectivity is essential for a satisfactory reconstruction. 

The regularized solution Bg is founded only on a suitable analysis of the real 
data, that aims at holding only the data whose information content is significant. 
This approach naturally agrees with the methodology of the experimental physical 
science. 

A moderate number of reasonable assumptions have been made in the construc- 
tion of the regularized solution Bg (see Theorem 13.511 . and, more important, the 
solution itself does not depend on unknown parameters. Even in the variational ap- 
proach, methods to reduce the dependence of the solution on free parameters have 
been widely investigated, and several practical strategies for choosing the regulariza- 
tion parameter a (see functional (|1.9|l ) have been proposed (see, for instance, |51 125j 
and references therein). Since the optimal parameter is impossible to determine be- 
cause the exact solution is not known, many of these strategies can provide estimates 
of the asymptotically optimal rate of convergence of the regularized solution to the 
real solution when the noise vanishes. 

The main difficulty of the method we have proposed regards the analysis of the 
correlation function. First, the correctness of the regularized solution depends on the 
capability of the correlation function to catch the information content of the data 
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and to exhibit it in an effective way. Second, usually quite large data samples, i.e., 
iV large, are necessary in order to limit sample fluctuations that could give rise to 
incorrect interpretation of the correlation function itself. 
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